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Abstract 

In the present work, we conduct large-scale orbital-free DFT calculations to study the en¬ 
ergetics of vacancy clustering in aluminum from electronic structure calculations. The sim¬ 
ulation domains considered in this study are as large as those containing a million atoms 
to accurately account for both the electronic structure and long-ranged elastic fields. Our 
results indicate that vacancy clustering is an energetically favorable mechanisms with pos¬ 
itive binding energies for a range of vacancy clusters considered in the present study. In 
particular, the 19 vacancy hexagonal cluster lying in {111} plane has a very large bind¬ 
ing energy with the relaxed atomic structure representative of a prismatic dislocation loop. 
This suggests that vacancy prismatic loops as small as those formed from 19 vacancies are 
stable, thus providing insights into the nucleation sizes of these defects in aluminum. 

Key words: Vacancy clustering, Prismatic dislocation, Electronic structure, Real space, 
Dislocation nucleation 


1 Introduction 


Prismatic dislocation loops play an important role in influencing the macroscopic me¬ 
chanical properties of materials, in particular ductility and fracture toughness in met¬ 
als (Gouldstone et al., 2000; Lubarta et al., 2004; Wirth, 2007). A large concentration of 
prismatic loops have been observed in quenched metals and in materials subject to large 
doses of radiation. It is widely believed that vacancy clustering is a precursor mechanism 
to the nucleation of prismatic loops in quenched metals ( Cuhlmann et al. , 1960; 2ot- 
terill & Segalla, 1963). In irradiated materials, experimental studies have shown a direct 
correlation between the dose of irradiation, the population of prismatic dislocation loops, 
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and the loss of ductility and fracture toughness (cf. e.g. fames & Mazey (196< ); Eyre & 
Bartlett (1965); Masters (1965); Trinkaus et. al. (1997); Singh et. al. (1997); Rice & Zin- 
kle (1998); Zinkle et. al (201 )). Numerous atomistic simulations have been conducted 
using empirical interatomic potentials to study the nucleation and evolution of defects in 
irradiated materials (cf. e.g. Bacon et al. (1993); Robinson (1994); Bacon & de la Rubia 
(1994); Ackland et al. (1997); Soneda & de la Rubia (1998); Trachenko et. al (2001); 
Wirth et al. (2000); Han et. al (2003); Marian et al. (2002a, b); Caturla et al. (2006)). These 
studies revealed displacement cascades nucleating a large concentration of vacancies and 
self-interstitals, which subsequently result in the formation of prismatic dislocation loops 
(vacancy loops and self-interstital clusters), among other defects, through the coalescence 
of vacancies into vacancy clusters and self-interstitials into self-interstitial clusters. 

While the displacement cascade simulations have elucidated the overall mechanism of the 
formation of prismatic loops, they do not provide insights into the energetics of formation 
of these defects which dictates the nucleation and evolution of these defects. Further, the 
nucleation of prismatic loops is a very rapid process that is hard to capture experimen¬ 
tally, especially given the high mobility of nanometer-sized prismatic loops (Arakawa et 
al., 2007; Matsukawa & Zinkle, 2007). To this end, atomistic calculations have been em¬ 
ployed to study the energetics of various sizes of vacancy and self-interstitial clusters (cf. 
e.g. Wirth et al. (1997, 2000); Marian et al. (2002a); Morishita et al. (2003)). While these 
atomistic studies have provided many important qualitative insights, these efforts have 
used empirical potentials to model interactions between various atoms in materials. The 
results of these simulations are significantly influenced by the choice of empirical poten¬ 
tials used and the physical quantities to which the parameters of these empirical potentials 
are fitted. It is hard to ascertain the accuracy of these empirical potentials to describe the 
details of the defect core, which involves situations with making and breaking of chemical 
bonds governed by quantum-mechanical interactions. Thus, it is desirable to conduct den¬ 
sity functional theory (DFT) calculations to study the energetics of the formation of these 
defects. The most widely used implementations of density functional theory are based on 
Fourier space formulations with a plane-wave discretization. While these Fourier space 
implementations have provided tremendous insights into the bulk properties of a wide 
range of materials, they are often restricted to simulation domains containing a few hun¬ 
dred atoms with periodic geometries. This poses a severe restriction in the electronic 
structure study of defects which require non-periodic geometries and larger simulation 
domains to accurately account for both the quantum-mechanical interactions at the defect 
core as well as the long-ranged elastic fields. 

In the present work we employ a real-space formulation of orbital-free density functional 
theory and the quasi-continuum reduction technique to conduct electronic structure calcu¬ 
lations on multi-million atom systems. In particular, we employ the Wang Govind Carter 
(WGC) orbital-free kinetic energy functional ( Vang et al , )) which has been shown 

to be accurate for aluminum for a wide range of properties ( darling & Carter, 20( ; Das 
et al , 5). We begin our study by computing the binding energies of divacancies in 

aluminum. A cell-size study of the binding energy of divacancies has shown significant 
cell-size dependence up to 2,000 atoms, underscoring the need for large-scale electronic 
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structure calculations to accurately study the energetics of defects in materials. The com¬ 
puted binding energies of both (100) and (110) divacancies are positive indicating the 
tendency of vacancies to attract. To understand the energetics of vacancy clustering, we 
next consider quad-vacancy clusters formed from divacancies. The computed binding en¬ 
ergies of all quad-vacancy clusters considered in this study are positive. Among the planar 
quad-vacancy clusters, those lying in the {111} plane had the highest binding energy. In 
order to investigate the energetics of vacancy clusters in the {111} plane, which is also 
one of the habit planes for vacancy prismatic loops, we consider hexagonal vacancy clus¬ 
ters of various sizes on this plane. While the 7 vacancy hexagonal cluster has positive 
binding energy, this is only marginally greater than the quad-vacancy cluster in the {111} 
plane. However, the 19 vacancy cluster has a very large relaxed binding energy and the 
atomic structure closely resembles a collapsed prismatic dislocation loop. This study sug¬ 
gests that vacancy clusters containing as small as 19 vacancies can collapse to form stable 
prismatic loops, thus providing insights into the nucleation size of these defects. 

The remainder of this paper is organized as follows. Section 2 presents a brief overview 
of the local real-space formulation of orbital-free DFT and the quasi continuum reduction 
technique employed in this work that has enabled consideration of multi-million atom 
computational domains. Section 3 presents our electronic structure study of vacancy clus¬ 
tering and nucleation of prismatic dislocation loops, along with a discussion of the new 
findings and their implications. We finally conclude with an outlook in Section 4. 


2 Overview of Methodology 


In this section, for the sake of completeness and to keep the discussion self-contained, we 
provide an overview of the local real-space formulation of orbital-free DFT, finite-element 
discretization, and the coarse-graining technique—quasi-continuum orbital-free DFT— 
employed in this work. We refer to Radhakrishnan & Gavini (2010); Motamarri et al. 
( 2) for a comprehensive discussion on the local real-space formulation, and Tavini et 

al. (2007a); Radhakrishnan & Gavini ( 010) for details on the coarse-graining technique. 
The main ideas are discussed below. 


2.1 Local real-space formulation of orbital-free DFT 


The ground-state energy in orbital-free density functional theory (Parr & Yarn , 200' ) is 
given by 


E(u, R) — T s ( u ) + E xc ( u ) + Eh ( u ) + E ext ( u , R) + E zz (R) , (1) 

where, u denotes the square-root of electron-density and R = {R 1; R 2 ,. .., Rm} de¬ 
notes the vector collecting the positions of atoms. In equation (1), T s denotes the kinetic 
energy of non-interacting electrons, which is explicitly modeled in orbital-free DFT; E xc 
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denotes the exchange-correlation energy, which includes all the quantum-mechanical in¬ 
teractions between electrons; E H denotes the Hartree energy or the classical electrostatic 
interactions between the electrons; E ext denotes the classical electrostatic interaction en¬ 
ergy between the electrons and nuclei; and E zz denotes the nuclear-nuclear electrostatic 
repulsion energy. 

The density-dependent Wang-Govind-Carter (WGC) ki netic energy functional ( Vang et 
a , 999) is the most widely used model for T s in solid-state orbital-free DFT calculations 
on materials systems whose electronic structure is close to a free-electron gas, and is 
employed in this work. Numerical studies have indicated that this is a transferable model 
for Al, Mg and Al-Mg materials systems with accuracies commensurate with Kohn-Sham 
DFT calculations on a wide range of materials properties (Carling & Carter, 2003; Ho et 
al., 2007; Das et al , 2015). The functional form of the WGC ki netic energy functional is 
given by 

T s (■ u) =C F J u 10/i dr + l,J | Vu (r)| 2 dr + T^(u) , (2) 

where 

T^ ,l3 (u) = C F J J u 2a (r)K(\r — r'|; w(r), u{r'))u 2 ^{r')drdr' . (3) 

In the above equation, the first term is the Thomas-Fermi energy with Cp = ^(37r 2 ) 2 / 3 , 
the second term is the von-Weizsacker correction (Pan - & Yang, 2003), and the last term, 
Tj" , denotes the density-dependent kernel kinetic energy functional. The kernel K and 
parameters a and /3 are chosen such that the linear response of uniform electron gas 
matches the theoretically known Lindhard response (Finnis, ). In particular, in the 
WGC functional, these parameters are chosen to be {a, (3} = {5/6+a/5/6,5/6-\/5/6}, 
and the density-dependent kernel is expanded as a Taylor series about the bulk average 
electron density, resulting in a series of density independent kernels (Wang et al., 1999). 


Widely used models for the exchange-correlation energy, especially for solid-state cal¬ 
culations of ground-state properties, include the local density approximation (LDA) and 
the generalized gradient approximation (GGA) (cf. e.g. Martin (2011)), and have been 
demonstrated to be transferable models for a range of materials systems and materials 
properties. In the present work, we will employ the LDA exchange-correlation functional 
given by 

E xc (u) = J e xc {u)u 2 {r)dr , (4) 

where e xc = e x + e c is the exchange and correlation energy per electron given by 


£ x (u) 


3 

4 




(5) 


j - jL - 

£ c (u) = < d+0iyf{rs)+02rs) * ~ ( 6 ) 

[ A log r s + B + C r s log r s + D r s r s < 1, 

where r s = (j^-) 1 / 3 . The values of constants used in this study are those of an unpolar¬ 
ized medium, and are given by = -0.1471, /3 lu = 1.1581, /3 2u = 0.3446, A u = 0.0311, 
B u = -0.048, C u = 0.0014, D u = -0.0108. 
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The remaining terms in equation (1) constitute the classical electrostatic energy between 
electrons and the nuclei, and are given by 



where Zj and V^(r, R j) denote the valence charge and pseudopotential corresponding 
to atom J located at R j, respectively. 

In the energy functional (1), all the terms are local, excepting the extended interactions in 
electrostatic energy and the kernel energy. A local variational reformulation of these ex¬ 
tended interactions in real-space has been developed in prior works ( iavini et al , 2007c; 
Radhakrishnan & Gavini, 201 0), which has enabled the consideration of general boundary 
conditions for studies on energetics of isolated defects (iyer et a] , 2015; Das et al., 1015) 
and is also a crucial step in developing coarse-graining schemes for electronic structure 
calculations. The extended interactions in electrostatics are governed by the j^ry ker¬ 
nel, which is the Green’s function of the Laplace operator. Thus, the total electrostatic 
interaction energy can be reformulated as the following local variational problem: 

E H +E ext +E zz = - inf j- 1 - j |V0(r)| 2 dr - J (u 2 ( r) + ^ 6/(r,Rj))0(r)dr j-E^/ • 

( 10 ) 

In the above, bj( r, R/) denotes the nuclear charge distribution corresponding to the ionic 
pseudopotential of the I th nucleus, </> denotes the electrostatic potential corresponding to 
the total charge distribution comprising of the electrons and the nuclei, and y is an ap¬ 
propriate function space. E set f denotes the self energy of the nuclear charge distributions 
which is computed by taking recourse to the Poisson equation associated with nuclear 
charge distribution (cf. Totamarri et al ( , 13)). 

Choly & Kaxiras (2002) proposed an approach to develop a local real-space formula¬ 
tion for the extended interactions in the kernel kinetic energy functional. In particular, 
they demonstrated that the series of density independent kernels obtained from the Tay¬ 
lor series expansion of the WGC density dependent kernel can each be approximated 
in the Fourier-space using a sum of partial fractions. Using this approximation, the ex¬ 
tended interactions in the kernel kinetic energy functional can be reformulated in terms 
of the solutions of a series of Helmholtz equations. If K (|r — r'|) denotes a density 
independent kernel, and is approximated in the Fourier-space using the approximation 
K (q) « HjLi , the kernel energy corresponding to the density independent kernel 

can be expressed as the following variational problem (Radhakrishnan & Gavini, 2010; 
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Motamarri et al., 2012): 


J J u 2ot (r)K(\r — r'\)u 2 P(r')drdr' 


inf sup L(u,u Q ,ojp), 

id a. £-2 id ft G-2 


(ID 


where 



In the above, ui a = {u > ai , co a2 ,, u am } and up = {c jp 1 , up 2 ,..., up m } denote the vector 
of potential fields, and Z denotes a suitable function space. We refer by ‘kernel poten¬ 
tials’ the auxiliary potential fields, u aj and up j for j = 1... m, introduced in the local 
reformulation of the extended interactions in the kernel energy. 


The problem of computing the ground-state energy and ground-state electronic structure 
of an N e electron system in orbital-free DFT for a fixed position of atoms can be expressed 
as a local variational problem in real-space using the reformulations in equations (10) and 
(11), and is given by 

inf sup inf sup L(u, </>, u a , up) subject to: / u 2 (r)dr = N e , (13) 

U£X 4>&y Ula&z U>p£Z J 

where L denotes the Lagrangian resulting from the local reformulations in equations 
(10) and (11). The function spaces, X, y and Z are appropriately chosen based on the 
boundary conditions dictated by the problem. The well-posedness of this local variational 
saddle-point orbital-free DFT problem has been established using the direct method in 
calculus of variations for certain models of orbital-free DFT kinetic energy functionals, 
and we refer to iavini et al. (2007c) for further details on the mathematical properties of 
the formulation. 


2.2 Finite-element discretization 


A finite-element basis presents a natural basis set to discretize the local real-space vari¬ 
ational formulation of orbital-free DFT discussed in section 2.1, and is employed in this 
work. The finite-element discretization of the orbital-free DFT problem presents many 
advantages over the widely used plane-wave discretization of orbital-free DFT calcu¬ 
lations ( ung et al., 2010). It allows for the consideration of complex geometries and 
boundary conditions that are not accessible through Fourier-space formulations of orbital- 
free DFT employing plane-wave discretization. This freedom from periodic boundary 
conditions is significant in the study of the energetics of defects in materials, which often 
break periodicity of perfect materials. To elucidate, the geometry of an isolated dislocation 
is not compatible with periodic boundary conditions which has limited our capabilities in 
studying the energetics of such defects. The recent developments in real-space formula¬ 
tion of orbital-free DFT and the finite-element discretization have enabled some of the 
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first studies on the energetics of an isolated edge dislocation in aluminum, providing key 
insights into the size of a dislocation-core or core-energetics from electronic structure 
calculations ( yer et al., 2015). While Fourier space implementations are computationally 
more efficient than real-space implementations using finite-element basis, recent numeri¬ 
cal efforts which use higher-order finite-element discretizations ( lotamarri et al., 2) 
have been shown to bridge this gap. Further, the good scalability of the finite-element 
discretization on parallel computing platforms allows the use of high performance com¬ 
puting resources to consider complex problems that are not accessible otherwise. Finally, 
the finite-element discretization is amenable to unstructured coarse-graining, which is the 
main idea exploited in the quasi-continuum coarse-graining technique, discussed subse¬ 
quently, which has enabled electronic structure calculations at macroscopic scales using 
orbital-free DFT. 


2.3 Quasi-continuum orbital-free density functional theory 


In this section, we present the key ideas in quasi-continuum orbital-free DFT (QC-OFDFT)— 
a seamless coarse-graining approach that enables us to perform multi-million atom orbital- 
free DFT simulations for electronic structure studies on defects at physically realistic 
concentrations. QC-OFDFT uses orbital-free DFT as the sole input physics and exploits 
the local real-space formulation in conjunction with the basis set adaptivity of the finite- 
element basis to achieve orbital-free DFT electronic structure calculations at macroscopic 
scales. The quasi-continuum reduction technique for orbital-free DFT was first demon¬ 
strated for local kinetic energy functionals ( lavim et al., 2007a) and later extended to 
the non-local WGC kinetic energy functionals employed in this work (Radhakrishnan & 
Gavini, 2010). 

The idea of quasi-continuum reduction was originally proposed in the context of inter¬ 
atomic potentials as a seamless scheme bridging the atomistic and continuum length- 
scales for crystalline materials (Tadmor et al., 1996; Knap & Ortiz , 2001) in order to 
simultaneously account for the atomistic interactions at the defect-core as well as the 
long-ranged elastic fields accompanying the defect. The quasi-continuum reduction idea 
is based on kinematic constraints introduced of the positions of atoms, thus reducing the 
number of degrees of freedom in the variational problem of computing the ground-state 
of a materials system. The kinematic constraints are introduced via a finite-element tri¬ 
angulation of a subset of atoms in the computational domain, which are referred to as 
the representative atoms or rep atoms. The rep atoms are chosen such that full atomistic 
resolution is provided in regions of interest, such as the defect-core with large atomistic 
displacements, and coarse-graining elsewhere. The kinematic constraints on atomistic dis¬ 
placements introduced through the finite-element triangulation of the rep atoms provides 
an excellent subspace for solving the variational problem of computing the ground-state 
energy of crystalline materials systems with defects. 

The quasi-continuum reduction for electronic structure calculations presents an additional 
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challenge of coarse-graining the electronic fields that exhibit oscillations on a subatomic 
length-scale. In QC-OFDFT formalism, the atomic displacements are treated in a simi¬ 
lar manner as the quasi-continuum approach for interatomic potentials where kinematic 
constraints are introduced on atomistic displacements using a finite-element triangula¬ 
tion of rep atoms. The electronic fields comprising of the square-root electron-density, 
electrostatic potential, and kernel potentials are decomposed into predictor and corrector 
fields: 


U =Uq + U c , 

0 =00 + 0C ) 

UJ a =UJ a o + UJac ) 

=UJ P0 + W /3 C > 


(14) 


where (u 0 , 0 O , u a0 , up 0 ) denote the predictor electronic-fields that are computed using 
periodic unit-cell calculations undergoing the Cauchy-Born deformation of the under¬ 
lying lattice. The predictor fields provide a good representation of the electronic fields 
to the leading order in regions of smooth deformations ( 31anc et al., 2002), such as re¬ 
gions away from the defect-core that are governed by elastic interactions. However, the 
electronic structure deviates significantly from that of the predictor fields in regions of 
rapidly varying deformations, such as in the defect-core. These deviations are captured 
through the corrector fields (u c , 0 C , Lu ac , u>p c ), and are computed from the local variational 
saddle-point problem discussed in section 2.1. As the predictor fields are a good repre¬ 
sentation of the electronic fields away from the defect-core, the corrector fields are signif¬ 
icant only near the defect-core. Thus, the corrector fields can be effectively represented 
by a finite-element triangulation which has sufficient resolution near the defect-core, but 
coarse-grains far away. The numerical implementation of the QC-OFDFT method is con¬ 
ducted using a hierarchy of finite-element triangulations (cf. figure 1): (i) an atomic-mesh 
representing atomic displacements, which has full atomistic resolution in the vicinity of 
the defect-core and coarse-grains far away; (ii) electronic-mesh representing the correc¬ 
tor electronic fields, which has subatomic resolution in the vicinity of the defect-core and 
coarse-grains away becoming superatomic; (iii) auxiliary unit-cell meshes used to repre¬ 
sent the predictor electronic fields. Often, for the sake of convenience, the electronic-mesh 
is chosen to be a subgrid of the atomic mesh. Using the decomposition in equation (14), 
the problem of computing the ground-state electronic structure reduces to a saddle-point 
variational problem in the coarse-grained independent variables comprising of the correc¬ 
tor fields on the electronic-mesh. We refer to davim et al (2007a) for further details and 
a comprehensive discussion of the method. 


Numerical studies have shown that the QC-OFDFT method converges rapidly with re¬ 
spect to coarse-graining, where a few thousand rep atoms have been sufficient to represent 
computational domains nominally containing millions of atoms for studying energetics 
of point defects. The effectiveness of QC-OFDFT has made possible electronic structure 
calculations at macroscopic scales using orbital-free DFT for studying the energetics of 
defects at realistic concentrations. Further, this provides a seamless scheme to account 
for both the quantum-mechanical effects at the defect-core and the long-ranged elastic 
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Fig. 1. (a) Atomic mesh used to represent positions of nuclei, (b) Auxiliary unit-cell meshes used 
to represent the predictor fields, (c) Electronic-mesh used to represent the corrector fields, which 
has subatomic resolution in the defect-core and coarsens away from the defect-core becoming 
superatomic. 

fields using an electronic structure theory (orbital-free DFT) as the sole input physics. 
Prior investigations have used QC-OFDFT to study the energetics of vacancies and dislo¬ 
cations in aluminum, and have provided important physical insights ( lav ini et al , 20071 ; 

Gavini, 2008; Radhakrishnan & Gavini, 2010). We refer to Iyer & Gavini (2011); Gavini 
& Liu (201 1) for numerical and mathematical analysis of the quasi-continuum reduction 
of orbital-free DFT, and field theories, in general. 


3 Results and Discussion 


We now present the results of our study on the energetics of vacancy clustering in alu¬ 
minum and nucleation of vacancy prismatic dislocation loops using large-scale orbital- 
free DFT electronic structure calculations. In this study, we employ the Wang-Govind- 
Carter (WGC) kinetic energy functional (Wang et al., 1999) with first-order Taylor ex¬ 
pansion of the density dependent kernel (cf. Wang et al. (1999)), a local density approx¬ 
imation (LDA) for the exchange-correlation energy ( ’erdew & Zunger, 198 ), and bulk 
local pseudopotential for aluminum (BLPS) ( hou et. al , 2( 4). The WGC ki netic energy 
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functional and the BLPS pseudopotential have been shown to be accurate and transfer¬ 
able for a range of material properties of Al, Mg and Al-Mg materials systems (Das et ah, 
2015). Numerical parameters such as finite-element discretization, coarsening of finite- 
element triangulation in the quasi-continuum method, numerical quadratures, and stop¬ 
ping tolerances on iterative non-linear and linear solvers are chosen such that the errors in 
the computed formation energies do not exceed 0.01 eV. In all the simulations reported 
in this work, we employ homogeneous Dirichlet boundary conditions on the corrector 
electronic fields, which correspond to the perturbations in the electronic fields arising 
from the defect vanishing on the boundary of the simulation domain with the electronic 
structure beyond the computational domain corresponding to that of the bulk. We refer 
to these boundary conditions on electronic fields as bulk Dirichlet boundary conditions 
which simulate an isolated defect in a bulk crystalline material. We hold the positions of 
atoms fixed on the boundary of the simulation cell, while relaxing the positions of the 
interior atoms. 


3.1 Divaccmcies 


In order to understand cell-size effects and establish the simulation domain sizes needed 
to accurately understand the energetics of vacancy interactions, we first conduct a cell-size 
study on divacancy binding energies in aluminum. The binding energy of an //-vacancy 
cluster is computed as: 

E'ZZ d = nE< v -El nv , (15) 

where E£ v is the formation energy of the //—vacancy cluster and E-f, is the formation 
energy of the monovacancy. The formation energy (at constant volume) of an n —vacancy 
cluster is given by (Finnis, 2003) 

El = E{N-n,n, t lj^v'}- t lj^E(N, 0, V), (16) 

where E (n — n. n. ;V Y n Vj denotes the energy of the system with N lattice sites oc¬ 
cupied by the n-vacancy cluster and N — n atoms with the total volume of the system 
being V. In the above, E (N. 0, V) denotes the energy of an iV-atom perfect crys¬ 
tal occupying a volume V. We compute the formation energies and binding energies at 
the equilibrium volume of a perfect crystal as we are primarily interested in the energet¬ 
ics of vacancy interactions in stress-free solids. We note that recent studies indicate that 
macroscopic deformations and macroscopic stresses influence the energetics of defects in 
a significant way (Gavini, 2008; Iyer et al., 2014, 2015), but this is not the focus of the 
present study and will be a topic for future investigations. 

Figure 2 shows the computed binding energies of divacancies along (100) and (110) crys¬ 
tallographic directions for cell-sizes ranging from 32 lattice sites to over 100,000 lattice 
sites. The results suggest a significant cell-size dependence on the computed binding en¬ 
ergies. In particular, this study indicates that cell-sizes containing about 2, 000 nominal 
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Fig. 2. Cell-size study of binding energies of divacancies along (100) and (110) directions in 
aluminum. 

number of atoms are required to obtain convergence in the binding energies of these sim¬ 
ple defects. These results underscore the long-ranged nature of the electronic fields and 
elastic displacement fields in the presence of defects, and emphasize the need for coarse- 
graining techniques such as QC-OFDFT for studying the energetics of defects. We note 
that the observed cell-size effects are in agreement with prior QC-OFDFT cell-size studies 
on mono- and di-vacancies formation energies in aluminum ( ravini et al., 2007a; Rad¬ 
hakrishnan & Gavini, 2010). The converged (100) and (110) divacancy relaxed binding 
energies are computed to be 0.045 eV and 0.067 eV, respectively. The positive binding 
energies of the divacancies suggest that vacancies tend to attract and form a stable diva¬ 
cancy complex as opposed to remaining as separated monovacancies. 


3.2 Quad-vacancies 


Having established the stability of divacancies, we next proceed to further investigate the 
vacancy clustering mechanism by computing the binding energies of quad-vacancy clus¬ 
ters formed from divacancies. As the number of possible quad-vacancy clusters is very 
large, we restrict our study to those quad-vacancy configurations where each vacancy 
has two other vacancies as the nearest or second nearest neighbors. This results in nine 
configurations, six of which are planar quad-vacancy clusters and three are non-planar 
configurations. Table 1 lists these configurations and the computed relaxed binding en- 
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ergies. In all these simulations, informed from the cell-size study of divacancies, we use 
cell-sizes containing 16,384 nominal number of atoms to ensure convergence with respect 
to cell-size. 

Table 1 

Binding energy of quad-vacancies in aluminum 


Structure 

Vacancy positions 

Binding energy (eV) 

Planar {100} 

(0,0,0), (a/2,a/2,0), (a,0,0), (a/2,-a/2,0) 

0.08 

Planar {100} 

(0,0,0), (a/2,a/2,0), (a,0,0), (3a/2,a/2,0) 

0.1 

Planar {100} 

(0,0,0), (a/2,a/2,0), (a,0,0), (a,a,0) 

0.1 

Planar {100} 

(0,0,0), (a,0,0), (0,a,0), (a,a,0) 

0.22 

Planar {110} 

(0,0,0), (0,a/2,a/2), (a,0,0), (a,a/2,a/2) 

0.3 

Planar {111} 

(0,0,0), (0,a/2,a/2), (a/2,a/2,0), (a/2,a,a/2) 

0.33 

Non Planar 

(0,0,0), (0,a/2,a/2), (a/2,0,a/2), (a/2,a/2,0) 

0.55 

Non Planar 

(0,0,0), (a,0,0), (a/2,a/2,0), (a/2,0,a/2) 

0.31 

Non Planar 

(0,0,0), (a,0,0), (a/2,a/2,0), (0,a/2,a/2) 

0.23 



Fig. 3. Electron-density contours of quad-vacancy cluster lying in {111} plane. 


The computed binding energies for all quad-vacancy clusters are positive suggesting an 
energetic preference to remain as quad-vacancy clusters as opposed to dissociation into 
monovacancies. However, the binding energies of the first three planar configurations of 
quad-vacancies lying in the {100} planes in table 1 are less than twice the binding energy 
of (110) divacancy, suggesting that these quad-vacancy clusters are likely to dissociate 
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Fig. 4. Electron-density contours of quad-vacancy cluster lying in {110} plane. 

into a pair of divacancies. On the other hand, the binding energies of other planar and 
non-planar quad-vacancy clusters are much larger and are stable with respect to dissocia¬ 
tion into divacancies. Among the planar quad-vacancy clusters, those lying in the {111} 
and {110} planes are the most stable. Figures 3 and 4 show the electron-density con¬ 
tours of quad vacancies lying in the {111} and {110} planes, respectively. Experimental 
studies suggest that these planes are the most likely habit-planes for vacancy clustering, 
eventually resulting in the formation of dislocation loops ( hlmann et al , 96'); Got¬ 
ten 11 & Segalla, 1963; Wu et al., 2007). The non-planar configurations are tetrahedral 
quad-vacancy clusters, with the largest binding energy corresponding to a stacking fault 
tetrahedra. Stacking fault tetrahedra formed from larger vacancy clusters play an impor¬ 
tant role in influencing plastic deformation in materials ( ■Ciritani, 1997; Kiritani et al., 
1999), and a study of the energetics of these defects will be a topic for future studies. In 
the present work, we will restrict ourselves to the study of larger planar vacancy clusters, 
in particular in the {111} plane, which is presented subsequently. 


3.3 Vacancy clusters in {111} and prismatic loops 


Experimental investigations indicate {111} as one of the potential habit planes for pris¬ 
matic dislocation loop nucleation from vacancy clustering ( ihlmann et al , ; Cot- 

terill & Segalla, 1963). It is hypothesized that when sufficiently large number of vacancies 
form a vacancy cluster, the planes of atoms above and below the cluster can collapse to 
form dislocation loops. In particular, studies indicate that dislocation loops with hexago¬ 
nal symmetry are commonly observed and are energetically favorable ( uhlmann et al. , 
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1960). However, the observed sizes of the loops are often 50 A diameter or larger, which 
are formed from vacancy clusters containing hundreds of vacancies. While loops of these 
sizes have been experimentally observed, the nucleation of loops is a very rapid process 
that is hard to observe experimentally. 



Fig. 5. Electron-density contours of unrelaxed hexagonal 7 vacancy cluster in {111} plane. 

In this work, we study the energetics of {111} hexagonal vacancy clusters to identify the 
nucleation size of vacancy prismatic loop in aluminum. As these are larger vacancy clus¬ 
ters, we employ a million atom simulation domain to accurately capture the electronic 
perturbations as well as the long-ranged elastic fields. We first begin with the study of the 
energetics of a 7-vacancy hexagonal cluster. Figure 5 shows the electron-density contours 
of an unrelaxed 7-vacancy hexagonal cluster in {111} plane. The relaxed binding energy 
of this vacancy cluster is computed to be 0.29 eV. The positive binding energy suggests 
that this vacancy cluster is stable with respect to dissociation into monovacanices. How¬ 
ever, the binding energy of this vacancy cluster is similar to the binding energy of quad- 
vacancy clusters on {111} planes. Thus, this vacancy cluster can potentially dissociate 
into smaller-sized vacancy clusters—for instance a quad-vacancy cluster, a divacancy and 
a monovacancy. Further, the relaxed configuration of this hexagonal vacancy cluster did 
not represent a collapsed vacancy loop. This is contrary to a prior study ( lavini et al., 
2007b), which showed a collapsed prismatic vacancy loop from a 7-vacancy hexagonal 
vacancy cluster. This discrepancy is primarily due to the Thomas-Fermi-von Weizsacker 
family of kinetic energy functionals employed in the prior study, which do not have the 
correct linear response of a free-electron gas. The WGC kinetic energy functional em¬ 
ployed in this work has the correct linear response of a free-electron gas, and benchmark 
studies have shown this kinetic energy functional to be in good agreement with Kohn- 
Sham DFT for a wide range of material properties in aluminum. 
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Fig. 6. Electron-density contours of unrelaxed hexagonal 19 vacancy cluster in {111} plane. 

We next study the energetics of the {111} hexagonal vacancy cluster containing 19 va¬ 
cancies. Figure 6 shows the electron-density contours of the unrelaxed hexagonal vacancy 
cluster. The relaxed binding energy of this hexagonal vacancy cluster is computed to be 
10.02 eV, which is a very large binding energy and significantly larger than that of the 7- 
vacancy cluster. The relaxed structure of this 19-vacancy hexagonal cluster resembles that 
of a collapsed prismatic loop. Figure 7 shows the relaxed atomic structure viewed along 
the (112) direction, which shows the collapse of planes to form a vacancy prismatic loop. 
Figure 8 shows the atoms in the simulation domain that do not have an face-centered-cubic 
(fee) coordination which depicts the core of the prismatic dislocation loop. The Burgers 
vector, based on the displacement fields of the collapsed loop, is calculated to be (| [111], 
Yg[112], [110]). The major component of the collapse of the planes surrounding the va¬ 

cancy loop is along (111) direction, which is perpendicular to the loop. In addition, the 
planes also slip along the (112) and (110) directions. The slip along the (112) direction is 
of particular importance. While the collapse along the (111) direction creates a prismatic 
dislocation loop, it also encloses a stacking fault which is energetically unfavorable. The 
slip along (112) eliminates the stacking fault, thus resulting in an unfaulted dislocation 
loop. The slip observed in the 19-vacancy loop is not sufficient to completely remove the 
stacking fault, but shows the tendency of the prismatic loop to unfault and we expect to 
observe complete unfaulting for larger-sized vacancy clusters. The present results suggest 
that vacancy clusters as small as 19 vacancies can collapse to form very stable prismatic 
dislocation loops, and represents the nucleation size of vacancy prismatic loops lying in 
the {111} planes. 
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Fig. 7. Atomic positions showing the collapse of vacancy cluster while viewing along (112). The 
atoms are color coded in the gray scale with depth along the viewing direction. 



Fig. 8. The core of the vacancy prismatic loop identified using coordination analysis. The depicted 
atoms are those in the simulation domain that are not in an fee coordination. 
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4 Conclusion 


We conducted large-scale electronic structure calculations by employing quasi-continuum 
orbital-free DFT to investigate the energetics of vacancy clustering in aluminum. In partic¬ 
ular, we used the WGC orbital-free kinetic energy functional in this study, whose accuracy 
has been ascertained for a wide range of properties in aluminum. The quasi-continuum 
reduction technique employed in this work has allowed electronic structure calculations 
using orbital-free DFT on multi-million atom simulation domains. Our cell-size study 
of binding energies of vacancies has shown that cell-sizes in excess of 2,000 atoms are 
required to obtain a converged value even for simple defects such as a divacancy. The 
binding energies of the divacancies and quad-vacancy clusters considered in this study 
are computed to be positive, suggesting the tendency of vacancies to coalesce to form va¬ 
cancy clusters. Among the planar quad-vacancy clusters the one on {111} had the highest 
binding energy, which is also a habit plane for the vacancy prismatic loops observed in ex¬ 
perimental investigations. Investigations on hexagonal vacancy clusters of varying sizes 
showed positive binding energies. However, the 19-vacancy hexagonal vacancy cluster 
had significantly higher binding energy compared to the 7-vacancy hexagonal cluster. An 
analysis of the relaxed atomic positions of the 19-vacancy cluster revealed a collapsed va¬ 
cancy prismatic dislocation loop. This study suggests that vacancy prismatic loops formed 
from vacancy clusters as small as 19 vacancies are stable and suggests that the nucleation 
sizes of these defects can be much smaller than the stable loops observed in experimental 
investigations that are often 50 A or larger in diameter. 

While this study has shown that vacancy clustering is an energetically feasible mechanism 
that can result in the nucleation of vacancy prismatic loops, many outstanding questions 
remain and are worthwhile topics for future investigations. An interesting and important 
question is relative stability of non-planar vacancy clusters of various sizes in comparison 
to planar vacancy clusters, which has important bearing on the formation of stacking fault 
tetrahedra. Another key question relates to the interaction of these defects with disloca¬ 
tions that influences the macroscopic mechanical properties of materials. Further, while 
orbital-free DFT provides important insights into the mechanisms and qualitative trends 
in the energetics, conducting similar studies using Kohn-Sham DFT will provide quanti¬ 
tatively accurate energetics. 
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